---
title: "16.Turnout"
---

# Load Packages
```{r}
library(tidyverse)
library(estimatr)
```

# Load Data
```{r}
load("Data/data_final.Rdata")
load("Data/gnr_panel.Rdata")
load("Data/turnout_list.Rdata")
```

# Set Covariates
```{r}
covs_list <- c("tx_Hi", "as.factor(DistN)", "ln_PopDens_70", "Cities_10km", "ln_Confs_All", "PrRecH_Adm3", "Isolados1950", "PrWheatAr_1972", "LandIneqRat", "TerrainSlopeIndex", "MedAltitude", "Communist_1975", "MinDistPCP_Map", "log_gnr_cumul_23_Z")
```

```{r}
elections <- c("AC75", "AR76", "AR79", "AR80", "AR83", "AR85", "AR87")
```

# Calculate First Differences in Turnout Over Time, Relative to 1975
```{r}

turnout_fd <- list()
for (i in 1:7){ # # for elections between 1976 and 1987
  T1 <- turnout_list[[1]]  # Get turnout for baseline year (1975)
  colnames(T1)[2:8] <- paste(colnames(T1)[2:8], "1", sep = "_")
  T2 <- turnout_list[[i]] # Get turnout for year of interest
  colnames(T2)[2:8] <- paste(colnames(T2)[2:8], "2", sep = "_")
  # First Difference
  turnout_fd[[i]] <- T2 %>% left_join(T1) %>% # Subtract 1975 turnout from turnout of year of interest
    mutate(Turnout = turnout_2 - turnout_1) %>%
    select(DICOFRE_Z, Turnout)
  is.na(  turnout_fd[[i]] ) <- do.call(cbind,lapply(  turnout_fd[[i]] , is.infinite))
}

```

#Cumulative
```{r}
  
baseline_list <- list()  
for (i in 2:7){ # i<-2

data_sample_g <- data_final %>% left_join(turnout_fd[[i]]) %>% # Join turnout first difference and cumulative GNR events up to election year to main data set
  left_join(gnr_panel %>% filter(elec_number == i)) 

output <- list()
for (g in 0:2){ # Subset for (3) groups (land reform intensity)
 
output[[g+1]] <- data.frame(election = elections[i], group = g,  model = colnames(turnout_fd[[i]])[-c(1)], tx_Hi = NA, se = NA, n = NA) %>% filter(model != "NA")

# Run turnout model
      m.1 <- lm_robust(as.formula(paste(output[[g+1]]$model, paste(c(covs_list), collapse= " + "), sep = " ~ ")),
                data = data_sample_g %>% filter(exp_cat2 == g), weight = weights, cluster = subclass)
      
           output[[g+1]][1, "tx_Hi"] <- m.1$coefficients["tx_Hi"]
           output[[g+1]][1, "se"] <- m.1$std.error["tx_Hi"]
           output[[g+1]][1, "n"] <- length(m.1$fitted.values)

}

output_list <- do.call(rbind, output)
baseline_list[[i]] <- output_list

}

# Plot Results

baseline_list_all <- do.call(rbind, baseline_list) 

years <- data.frame(election = c("AR76", "AR79", "AR80", "AR83", "AR85", "AR87"), 
           year = as.factor(c(1976, 1979, 1980, 1983, 1985, 1987)))
baseline_list_all <- baseline_list_all %>% left_join(years)

ggplot(baseline_list_all, aes(x = year, col = as.factor(group))) +
  geom_point(aes(y = tx_Hi, shape = as.factor(group)), position = position_dodge(width = 1)) +
  geom_linerange(aes(ymin=tx_Hi - (qt(0.975, n)*se), ymax=tx_Hi + (qt(0.975, n)*se)), position = position_dodge(width = 1), size = 0.5) +
  geom_linerange(aes(ymin=tx_Hi - (qt(0.95, n)*se), ymax=tx_Hi + (qt(0.95, n)*se)), position = position_dodge(width = 1), size = 0.75) +
  geom_hline(aes(yintercept = 0), col = "red") +
  facet_wrap(~model, scales = "free_y") +
  xlab("Election Year") + ylab("Treatment Effect") +
  scale_color_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                     values = c("#252525", "#636363", "#969696"))+
  scale_shape_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                       values = c(16, 17, 15)) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 90))

ggsave(plot = last_plot(),
         "Figures/Afg6a_RobustTurnout.pdf",
             width = 9, height = 5)
  
```

